//__gnu_debug::set
#include "pch.h"

#define NDEBUG = 1
#define BOOST_DISABLE_ASSERTS
//#define DEBUG_STL = 1

//
//  Standard headers
// 
#include <iostream>
	using std::cerr;
	using std::cout;
	using std::ofstream;
	using std::istream;
#include <utility>
	using std::make_pair;

#include <set>
#include <vector>

#include <string>
	using std::string;
#include <numeric>
#include <algorithm>
	//using std::fill_n;
#include <stdexcept>
	using std::runtime_error;
#include <unordered_map>
	using std::unordered_map;
	typedef unordered_map<unsigned, unsigned> map_uu;
#include <tuple>
	using std::tuple;
	using std::tie;

#include <open_stream_name_with_substitutions.h>


//
//  My common headers
// 
#include <progress_indicator.h>
#include <commify.h>
#include "portable_timer.h"
#include "intervals_per_contig.h"
	using intervals_per_contig::t_intervals_per_contig;
	using intervals_per_contig::t_interval_search_cache;
	using intervals_per_contig::t_intervals;
	using intervals_per_contig::ci_map_str_intervals;



//
//  Boost headers
// 
#include <boost/regex.hpp>
	using boost::regex;
#include <boost/algorithm/string/replace.hpp>
	using boost::replace_all_copy;
#include <boost/functional/hash.hpp>
	using boost::hash_combine;
#include <boost/iterator/counting_iterator.hpp>
	using boost::make_counting_iterator;
#include <boost/spirit/include/qi.hpp>
	namespace qi = boost::spirit::qi;
#include <boost/fusion/include/adapt_struct.hpp>
#include <boost/fusion/include/io.hpp>

#include <boost/filesystem.hpp>
	using boost::filesystem::path;
	using boost::filesystem::file_size;
	using boost::filesystem::is_regular_file;
	//using boost::filesystem::basename;

#include <boost/array.hpp>
	using boost::array;
#include <boost/foreach.hpp>

#include <boost/range/algorithm_ext/push_back.hpp>
    using boost::push_back;
#include <boost/range/algorithm.hpp>
	using boost::sort;
	using boost::reverse;
	using boost::for_each;
#include <boost/range/numeric.hpp>
	using boost::accumulate;

#include <boost/algorithm/string.hpp>
    using boost::join;
	using boost::algorithm::replace_all;
#include <boost/format.hpp>
	using boost::format;
#include <boost/iostreams/device/file.hpp>
#include <boost/iostreams/stream.hpp>
#include <boost/iostreams/device/mapped_file.hpp>
	namespace io = boost::iostreams;
#include <boost/assert.hpp> 

//#include <boost/spirit/include/support_multi_pass.hpp>

//#include <boost/spirit/include/karma.hpp>
//	namespace karma = boost::spirit::karma;
//	using boost::spirit::karma::uint_;
//	using boost::spirit::karma::lit;
//	using boost::spirit::karma::generate;

#include <boost/fusion/include/adapt_struct.hpp>
#include <boost/fusion/include/io.hpp>
#include <boost/spirit/include/phoenix_core.hpp>
namespace phoenix = boost::phoenix;
#include "qi_parse_error_handler.h"




#include "read_options.h"

#include <cassert>
#include "boost/multi_array.hpp"
	using boost::multi_array;
	using boost::const_multi_array_ref;
	using boost::extents;

#ifdef DEBUG_STL
	#include <debug/set>
	#include <debug/vector>
	using __gnu_debug::set;
	using __gnu_debug::vector;
#else
	using std::set;
	using std::vector;
#endif


void do_main(prog_options& opt);

//________________________________________________________________________________________
// 
//  	Main
//			boilerplate to catch errors
// 
//________________________________________________________________________________________
int main (int argc, char *argv[])
{
	try
	{
		// 
		//  Read options
		//  	return if user chose help or version
		// 
		prog_options opt;
		if (!read_options(opt, argc, argv))
			return 0;

		do_main(opt);

	}
	catch(const std::exception& e)
	{
		// indent error string
		string error_string = e.what();
		replace_all(error_string, "\n", "\n\t");
		std::cerr
				<< "\n"<< string(60, '=') << "\n"
				<< "\nError:\n\t" << error_string 
				<< "\n\n" << string(60, '=') << "\n\n";
		return -1;
	}
	catch (...)
	{
		std::cerr
				<< "\n"<< string(60, '=') << "\n"
				<< "\nError:\n\tUnexpected exception of unknown type\n\n"
				<< string(60, '=') << "\n\n";
		return -1;
	}
	return 0;
}



void throw_logic_error(format& f)
{
	throw std::logic_error(f.str());
}






//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//		Read founder alleles

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//!_______________________________________________________________________________________
//!     
//!		t_founder_contig
//!			alleles per contig
//!     
//!		\author lg (6/4/2011)
//!_______________________________________________________________________________________ 
class t_founder_contig
{
public:
	t_founder_contig():
		cnt_founder_strains (0),
		cnt_sites           (0),
		cnt_diallelic_sites (0)
		{}

	unsigned   					cnt_founder_strains;
	unsigned   					cnt_sites;
	unsigned   					cnt_diallelic_sites;

	string 						contig;
	vector<unsigned>			loci;
	vector<char>				ok;

								// snps are provided as a 2-d matrix but stored for reading as a vector
	const_multi_array_ref<char, 2>	get_diallellic_snps() const {assert (diallellic_snps_storage.size() == cnt_sites * cnt_founder_strains);
																return const_multi_array_ref<char, 2>(&diallellic_snps_storage[0],
																						extents[cnt_sites][cnt_founder_strains]);}
	vector<char>				diallellic_snps_storage;

	
	//							// alleles are provided as a 2-d matrix but stored for reading as a vector
	//const_multi_array_ref<string, 2>	get_all_alleles() const {assert (all_alleles_storage.size() == cnt_sites * cnt_founder_strains);
	//															return const_multi_array_ref<string, 2>(&all_alleles_storage[0],
	//																				extents[cnt_sites][cnt_founder_strains]);}
	//vector<string>				all_alleles_storage; 

								// lookup table converting from contig position to index into alleles
								// -1 if no or ignored allele
	vector<int>					map_pos_to_snp_index;
	const_multi_array_ref<char, 2>	get_sdp() const {assert (sdp_storage.size() == cnt_sites * cnt_founder_strains);
																return const_multi_array_ref<char, 2>(&sdp_storage[0],
																			extents[cnt_sites][cnt_founder_strains]);}
	vector<char>				sdp_storage;

	void						sanity_check() const
	{
		assert(loci.size() == cnt_sites);
		assert(ok.size() == cnt_sites);
		assert(diallellic_snps_storage.size() == cnt_sites * cnt_founder_strains);
		//assert(all_alleles_storage.size() == cnt_sites * cnt_founder_strains);
		assert(sdp_storage.size() == cnt_sites * cnt_founder_strains);
		assert(cnt_sites);
		assert(cnt_founder_strains);
	}
};



//!_______________________________________________________________________________________
//!     
//!		read_strain_names_from_header
//!     
//!		\author lg (6/4/2011)
//!_______________________________________________________________________________________
void read_strain_names_from_header(std::istream& input_stream, const string& file_name, vector<string>& strain_names)
{
	string line_str;
	getline(input_stream, line_str);
	auto beg = line_str.begin();
	qi::parse(beg, line_str.end(), +qi::graph % qi::omit[qi::space], strain_names);
	if (strain_names.size() <= 5)
		throw std::logic_error("No strain names in header of allele file [" + file_name + "]");
	if (strain_names[0] != "chr")
		throw std::logic_error("Ill formed header line does not begin with chr: [" + file_name + "]");
	
	// get rid of first five fields which are before the strain_names
	strain_names.erase(strain_names.begin(), strain_names.begin() + 5);
}




struct t_founder_site
{
	string			contig;
	unsigned		pseudo_pos;
	double   		pos;
	unsigned		cnt_alleles;
	unsigned		maf;
	vector<string>	alleles;
};

BOOST_FUSION_ADAPT_STRUCT(
t_founder_site,
(string        ,  contig)    
(unsigned      ,  pseudo_pos)
(double         ,  pos)       
(unsigned      ,  cnt_alleles)  
(unsigned      ,  maf)       
(std::vector<std::string>,  alleles)   
)

phoenix::function<errorHandlerT> const errorHandler_body = errorHandlerT("Founder alleles");

 
//	Parsing grammar 
// 
//  chr   pse.bp   bp     nalleles   maf   SPECIES_1   SPECIES_2   SPECIES_3    
//  1     11       11     2          18    T           T           T            
//  1     22       21.5   3          16    0           0           0            
template <typename Iterator>
struct t_founder_site_grammar
: qi::grammar<Iterator, t_founder_site()>
{
    t_founder_site_grammar(unsigned& curr_line_number_, string& error_msg)
    : t_founder_site_grammar::base_type(founder_site), curr_line_number(curr_line_number_)
    {
		using qi::graph;
		using qi::uint_;
		using qi::double_;
		using qi::char_;
		using qi::lit;
        using namespace qi::labels;
        using qi::eoi;
        using qi::eol;

		contig		%= -(lit("Chr") | lit("chr")) > +graph;
		pseudo_pos	%= uint_;
		pos			%= double_;
		cnt_alleles	%= uint_;
		maf			%= uint_;
		allele		%= +graph;
		whitespace  = qi::space;


		founder_site %=	 			 	contig 					// contig
						 > whitespace >	pseudo_pos				// pse.bp               
						 > whitespace >	pos       				// bp                   
						 > whitespace >	cnt_alleles				// nalleles             
						 > whitespace >	maf       				// maf                  
						 > whitespace > allele % whitespace		// strain_1 .. strain_N 
						 > eoi;

        // give names for debug
        contig      .name("Chromosome");
        pos         .name("Contig position");
		pseudo_pos  .name("Contig 'pseudo' position");
        cnt_alleles .name("Number of alleles");
        maf         .name("Major Allele Frequency");
        allele      .name("allele");

        qi::on_error<qi::fail>                               
        (founder_site, ::errorHandler_body(_1, _2, _3, _4, phoenix::ref(curr_line_number), phoenix::ref(error_msg), true));
    }
	unsigned& curr_line_number;
    qi::rule<Iterator, std::string()>     		contig;
    qi::rule<Iterator, unsigned()>      		pseudo_pos;
	qi::rule<Iterator, double()>      			pos;
	qi::rule<Iterator, unsigned()>      		cnt_alleles;
	qi::rule<Iterator, unsigned()>      		maf;
	qi::rule<Iterator, std::string()>       	allele;
	qi::rule<Iterator, void()>    				whitespace;
    qi::rule<Iterator, t_founder_site()>		founder_site;
};


//!_______________________________________________________________________________________
//!     
//!		read_allele_data
//!			read all the alleles for one chromosome
//!     
//!		\author lg (7/4/2011)
//!     
//!_______________________________________________________________________________________
void read_allele_data(	prog_options& 					opt, 
						const string& 					file_name, 
						vector<t_founder_contig>&		founder_contigs,
						vector<string>&					contig_names,
						vector<string>&					strain_names)
{
	opt.both() << "Parsing " + file_name + "\n";

	// 
	// guess number of rows for progress dots based on file sizes
	// 
	//path input_path(file_name);
	if (!is_regular_file(file_name))
		throw std::logic_error(file_name + " does not exist or is not a regular file.");

	// 120 bytes per line (135 empirically /data/www/18genomes/variants.tables/*)
	auto input_file_size = file_size(file_name);
	unsigned guess_rows = (input_file_size - 200) / 120 + 1;


	io::stream<io::mapped_file_source> input_stream;
	input_stream.open(file_name.c_str());
	if (!input_stream)
		throw std::logic_error("Can't open " + file_name);





	//	
	//	Add data for new contig
	//
	founder_contigs.push_back(t_founder_contig());
	t_founder_contig& ad = founder_contigs.back();


	//
	//	get strain names 
	// 
	// 
	vector<string>			strain_names_in_file;
	read_strain_names_from_header(input_stream, file_name, strain_names_in_file);
	// use these strain_names if not previously defined
	if (!strain_names.size())
	{
		strain_names.swap(strain_names_in_file);
		// log the list of strains
		opt.both() << "  " << strain_names.size() << " strains\n";
		opt.both(3)  << "  [" + join(strain_names, ", ") + "]\n";
	}
	else
		// check strain names are the same across chromosomes!!
		if (strain_names != strain_names_in_file)
			throw std::logic_error("Strain names in [" + file_name + "] = "
									"[" + join(strain_names_in_file, ", ") + "]"
									" are different from those previously parsed = "
									"[" + join(strain_names,         ", ") + "]");
	ad.cnt_founder_strains = strain_names.size();




	// 
	// ignore ambiguous alleles as well
	// 

	unsigned char_range = std::numeric_limits<char>::max() - std::numeric_limits<char>::min() + 1;
	vector<char> valid_nucleotides (char_range, 1);
	for_each(string("ACGTacgt"), [&](unsigned c){valid_nucleotides[c] = 0;});

	unsigned cnt_rows = 0;
	string error_msg;
	t_founder_site_grammar<string::const_iterator> grammar(cnt_rows, error_msg);

	unsigned cnt_biallelic = 0;
	unsigned cnt_all_snps  = 0;

	{
		// progress dots
		t_progress_indicator progress(opt.verbose(), guess_rows / 49 + 1, guess_rows);
		progress.use_timer();
	
		string line_str;
		while (getline(input_stream, line_str))
		{
			++progress;
			cnt_rows += 1;
	
			// 
			// 	Parse whitespace separated elements
			// 		Name unused fields just in case...
			// 
			t_founder_site	founder_site;
			auto beg = line_str.cbegin();
			auto end = line_str.cend();
			if (!qi::parse(beg, end, grammar, founder_site) || beg != end)
				throw_logic_error(format("Parsing error on Row %1$d.") % cnt_rows);
				//throw std::logic_error((format("Parsing error on Row %1$d.") % cnt_rows).str());
			//cerr	<< end -beg							<< ", " 
			//		<< founder_site.contig 				<< ", " 
			//		<< founder_site.pseudo_pos 			<< ", " 
			//		<< founder_site.pos 				<< ", "
			//		<< founder_site.cnt_alleles 		<< ", " 
			//		<< founder_site.maf 				<< "\n";
	
	
			// check right number of alleles
			if (founder_site.alleles.size() != ad.cnt_founder_strains)
				throw_logic_error(format("Expecting %1$d alleles on every row. "
											   "Row %2$d has %3$d values.") 
											% ad.cnt_founder_strains
											% cnt_rows
											% founder_site.alleles.size());
	
			ad.cnt_sites += 1;
	
			// chop off fractional position and save
			ad.loci.push_back(founder_site.pos);
	
	
			// add to all_alleles
			//push_back(ad.all_alleles_storage, founder_site.alleles);
			//ad.all_alleles_storage.insert(ad.all_alleles_storage.end(), founder_site.alleles.begin(), founder_site.alleles.end());
	
	
			// 
			// save contig name
			//
			if (!ad.contig.length())
			{
				contig_names.push_back(founder_site.contig);
				ad.contig = founder_site.contig;
			}
	
			//
			//	check if is diallelic ACGT SNP 
			// 
			char ok = 1;
			if (founder_site.cnt_alleles != 2)
				ok  = 0;
			else
			{
				++cnt_biallelic;
	
				// 
				//  All SNPs? (no indels):
				//  	 if total length = number of alleles (1 bp each)
				// 
				if (founder_site.alleles.size() != accumulate(founder_site.alleles, 0, [](int len, const string& a){return len + a.length();}))
					ok = 0;
				else
				{
					++cnt_all_snps;
	
					// 
					//  All ACGT? (no "0"):
					//  	 only ACGT = 0, i.e. sum of all alleles == 0 if all valid
					// 
					if (0 != accumulate(founder_site.alleles, 0, [&](int v, const string& a) {return v + valid_nucleotides[static_cast<unsigned char>(a[0])];}))
						ok = 0;
					else
						++ad.cnt_diallelic_sites;
				}
			}
				
	
			// save whether site is ok or not
			ad.ok.push_back(ok);
	
			//
			//	Diallelic valid SNPs 
			// 
			if (ok)
			{
				//  Copy diallelic SNPs (first letter of each allele) to diallellic_snps
				for_each(founder_site.alleles, [&](const string& a) {ad.diallellic_snps_storage.push_back(a[0]);});
	
				//  
				// Strain distribution pattern
				// 
				//  For diallelic SNPs must be either equal to first allele or not!
				char first_allele = founder_site.alleles[0][0];
				BOOST_FOREACH(const string& allele, founder_site.alleles) 
					ad.sdp_storage.push_back(allele[0] == first_allele ? 0: 1);
			}
	
			// 
			//  For alleles with possibly indels
			// 
			else
			{
				// hyphens for invalid sites
				ad.diallellic_snps_storage.insert( ad.diallellic_snps_storage.end(), ad.cnt_founder_strains, '-');
	
				//  
				// 	Strain distribution pattern
				//  	successive newly encountered alleles = incrementing indices
				// 
				
				
				// 
				// only change map_sdp_to_index if it wasn't there already
				//
				unordered_map<string, unsigned> map_allele_to_index;
				BOOST_FOREACH(const string& allele, founder_site.alleles) 
				{
					auto ret = map_allele_to_index.insert(make_pair(allele, map_allele_to_index.size()));
					ad.sdp_storage.push_back(ret.first->second);
				}
			}
	
		}
	}

	// 
	//  DEBUG
	// 
	//const_multi_array_ref<char, 2>	sdp			= ad.get_sdp();
	//for (unsigned xx = 0, xx_end = sdp.size(); xx < xx_end; ++xx)
	//{
	//	cout << ad.loci[xx];
	//	for (unsigned yy = 0, yy_end = sdp[0].size(); yy < yy_end; ++yy)
	//		cout << "\t" << static_cast<unsigned>(sdp[xx][yy]);
	//	cout << "\n";
	//}

	opt.both(2) << format("Chromosome %5s: %7s rows , #sites: biallelic = %7s / "
						 "SNPs = %7s / Valid SNPs = %7s\n") 
					% ad.contig
					% commify(cnt_rows)
					% commify(cnt_biallelic)
					% commify(cnt_all_snps)
					% commify(ad.cnt_diallelic_sites);
	ad.sanity_check();
}


//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//		Read short read SNP calls

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
//!_______________________________________________________________________________________
//!     
//!		t_seq_calls_per_contig
//!			short read calls per contig
//!     
//!		\author lg (7/4/2011)
//!_______________________________________________________________________________________ 
class t_seq_calls_per_contig
{
public:
		vector<unsigned>	pos;
		vector<char>		alleles;
		vector<unsigned>	indices_into_founder_contig;
		void				sanity_check() const
		{
			assert (pos.size() == alleles.size() );
			assert (pos.size() == indices_into_founder_contig.size());
		}

};



//!_______________________________________________________________________________________
//!     
//!		read_seq_calls
//!			read short read calls per strain
//!     
//!		\author lg (7/4/2011)
//!_______________________________________________________________________________________
void read_seq_calls(	prog_options& 							opt, 
						const string& 							file_name, 
						unsigned								cnt_founders,
						vector<map_uu>&							map_pos_to_allele_index_per_contig,
						vector<t_seq_calls_per_contig>&			seq_calls_per_contig,
						vector<string>&							contig_names)
{
	//BOOST_FOREACH(const string& s, contig_names) 
	//	opt.verbose() << s << "\n";
	

	// 
	// guess number of rows for progress dots based on file sizes
	// 
	//path input_path(file_name);
	if (!is_regular_file(file_name))
		throw std::logic_error(file_name + " does not exist or is not a regular file.");

	// 17 bytes per line for progress dots
	//auto input_file_size = file_size(file_name);
	//unsigned guess_rows = input_file_size / 17 + 1;

	io::stream<io::mapped_file_source> input_stream;
	input_stream.open(file_name.c_str());
	if (!input_stream)
		throw std::logic_error("Can't open " + file_name);





	//	
	//	Add data for each contig
	//
	seq_calls_per_contig.assign(contig_names.size(), t_seq_calls_per_contig());

	//
	//	map contig names to indices in alleles per contig
	// 
	unordered_map<string, unsigned> map_contig_to_indices;
	for (unsigned ii = 0; ii < contig_names.size(); ++ii)
		map_contig_to_indices[contig_names[ii]] = ii;


	// 
	//	Parsing grammar 
	// 
	using qi::graph;
	using qi::uint_;
	using qi::char_;
	using qi::lit;
	using qi::omit;
	//  Chr1	11	T	?
	auto contig = -(lit("Chr") | lit("chr")) > +graph;	// ignore leading "Chr"
	auto space = qi::omit[qi::space];
	auto grammar = 	contig			> space > 			// chr
					uint_			> space >   		// pos
					omit[char_] 	> space >   		// omit reference allele
					char_;								// alt allele

	// 
	// ignore ambiguous alleles
	// 
	unsigned char_range = std::numeric_limits<char>::max() - std::numeric_limits<char>::min() + 1;
	vector<char> invalid_nucleotides (char_range, 1);
	BOOST_FOREACH(unsigned c, string("ACGTacgt")) 
		invalid_nucleotides[c] = 0;


	vector<unsigned> cnt_alleles_in_founder_per_contig 	(contig_names.size(), 0);
	vector<unsigned> cnt_alleles_per_contig 			(contig_names.size(), 0);
	vector<unsigned> cnt_valid_alleles_per_contig 		(contig_names.size(), 0);

	{
		// Two rows of dots progress dots
		//t_progress_indicator progress(opt.verbose(), guess_rows / 48 + 1, guess_rows);
		//progress.use_timer();
	
		unsigned cnt_rows = 0;
		string line_str;
		while (getline(input_stream, line_str))
		{
			//++progress;
			cnt_rows += 1;
	
			//
			// 		Parse
			// 
			string 			contig;
			unsigned		pos;
			char			allele;
			auto beg = line_str.begin();
			auto end = line_str.end();
			qi::parse(beg, end, grammar, contig, pos, allele);
	
	
	
			if (!map_contig_to_indices.count(contig))
				throw std::logic_error("Chromosome " + contig + " not in founder allele files");
			unsigned contig_index = map_contig_to_indices[contig];

			++cnt_alleles_per_contig[contig_index];
	
			// ignore non ACGT
			if (invalid_nucleotides[static_cast<unsigned char>(allele)])
				continue;
			++cnt_valid_alleles_per_contig[contig_index];
	
			// Ignore alleles not in founder strains
			assert(contig_index < map_pos_to_allele_index_per_contig.size());
			auto ff =map_pos_to_allele_index_per_contig[contig_index].find(pos);
			if (ff == map_pos_to_allele_index_per_contig[contig_index].end())
				continue;
	
			seq_calls_per_contig[contig_index].pos.push_back(pos);
			seq_calls_per_contig[contig_index].alleles.push_back(allele);
			seq_calls_per_contig[contig_index].indices_into_founder_contig.push_back(ff->second);
			++cnt_alleles_in_founder_per_contig[contig_index];
		}
	}

	for (unsigned ii = 0; ii < contig_names.size(); ++ii)
	{
		opt.log() << format("Chromosome %s: %12s alleles / %12s valid alleles / %12s alleles in founder\n") 
						% contig_names[ii] 
						% cnt_alleles_per_contig[ii]
						% cnt_valid_alleles_per_contig[ii]
						% cnt_alleles_in_founder_per_contig[ii];
	}
}



//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//		Read intervals to exclude

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//!_______________________________________________________________________________________
//! 
//!		read_intervals
//!			read intervals in gff or tab delimited format
//! 
//!		\author lg (08/04/2011)
//!     
//!_______________________________________________________________________________________
void read_intervals(	prog_options& 										opt, 
						const string&										file_name, 
						t_intervals_per_contig& 							intervals_by_contig, 
						bool 												zero_based_coordinates, 
						intervals_per_contig::eINTERVAL_FILE_FORMAT 		file_format)
{
	opt.both() << "Reading intervals from " + file_name + "\n";
	io::stream<io::mapped_file_source> input_stream;
	input_stream.open(file_name.c_str());
	if (!input_stream)
		throw std::logic_error("Can't open " + file_name);

	intervals_per_contig::read_intervals(input_stream, intervals_by_contig, zero_based_coordinates, file_format, false, &opt.both(4));
}


//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//		Reconstruction

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//!_______________________________________________________________________________________
//! 
//!		t_founder_interval
//!			results of reconstruct mosaic
//! 
//!		\author lg (09/04/2011)
//!_______________________________________________________________________________________
struct t_founder_interval
{
	//unsigned			beg_index;
	//unsigned			back_index;
	unsigned			founder_contig_beg_index;
	unsigned			founder_contig_back_index;
	unsigned			beg_locus;
	unsigned			back_locus;
	unsigned			founder_strain_id1;
	unsigned			founder_strain_id2;
	unsigned			cnt_mismatches;
	unsigned			cnt_transitions;
	unsigned			cnt_informative;
	unsigned			cnt_uninformative;
	unsigned			cnt_snps;
	unsigned			len() const { return back_locus - beg_locus + 1;}
	double 				errors_per_snp() const { return cnt_snps > 0 ? cnt_mismatches / (double)cnt_snps : 0.0;}
	double 				errors_per_bp () const { return len()	 > 0 ? cnt_mismatches / (double)len() 	 : 0.0;}
	vector<unsigned>	mismatch_loci;

	// constructor
	t_founder_interval(	unsigned			founder_contig_beg_index_,
						unsigned			founder_contig_back_index_,
						unsigned			beg_locus_,
						unsigned			back_locus_,
						unsigned			founder_strain_id1_,
						unsigned			founder_strain_id2_,
						unsigned			cnt_mismatches_,
						unsigned			cnt_transitions_,
						unsigned			cnt_informative_,
						unsigned			cnt_uninformative_,
						unsigned			cnt_snps_):
								founder_contig_beg_index(founder_contig_beg_index_),          
								founder_contig_back_index(founder_contig_back_index_),         
								beg_locus (beg_locus_),
								back_locus(back_locus_),
								founder_strain_id1(founder_strain_id1_),
								founder_strain_id2(founder_strain_id2_),
								cnt_mismatches    (cnt_mismatches_),
								cnt_transitions   (cnt_transitions_),
								cnt_informative   (cnt_informative_),
								cnt_uninformative (cnt_uninformative_),
								cnt_snps(cnt_snps_)
								{
								}
	t_founder_interval(	unsigned	founder_contig_beg_index_,
						unsigned	founder_contig_back_index_,
						unsigned	beg_locus_,
						unsigned	back_locus_,
						unsigned	founder_strain_id_,
						unsigned	cnt_mismatches_,
						unsigned	cnt_snps_):
								founder_contig_beg_index(founder_contig_beg_index_),          
								founder_contig_back_index(founder_contig_back_index_),         
								beg_locus (beg_locus_),
								back_locus(back_locus_),
								founder_strain_id1(founder_strain_id_),
								founder_strain_id2(founder_strain_id_),
								cnt_mismatches    (cnt_mismatches_),
								cnt_transitions   (0),
								cnt_informative   (cnt_snps_ - cnt_mismatches_),
								cnt_uninformative (0),
								cnt_snps(cnt_snps_)
								{
								}
};


//!_______________________________________________________________________________________
//!     
//!		find_optimal_path_via_dp
//!			
//!     
//!		\author lg (4/11/2011)
//!     
//!     \param mat 					Matrix
//!     \param cnt_loci 
//!     \param matrix_size 
//!     \param boundaries 			Where founders change (including beginning and end)
//!     \param founder_ids 			Successive list of founders
//!_______________________________________________________________________________________
void find_optimal_path_via_dp(	const multi_array<char,2>& 		mat,
								unsigned						cnt_loci,
								unsigned						matrix_size,
								vector<unsigned>&				boundaries,
								vector<unsigned>&				founder_ids,
								const multi_array<unsigned,2>& 	penalty_matrix)
{
	assert(cnt_loci);
	assert(matrix_size);
	multi_array<int,2> path(extents[cnt_loci][matrix_size]);
	//multi_array<int,2> 	opt (extents[cnt_loci][matrix_size]);
	vector<int>			opt_curr(matrix_size, 0);
	vector<int>			opt_prev(matrix_size, 0);
	//
	//	initialise first column to matrix
	//
	for (unsigned jj = 0; jj < matrix_size; jj++)
	{
		opt_prev[jj] = mat[0][jj];
		path[0][jj] = jj;
	}

	// given previous column in opt, what is the best path to cells in opt
	for (unsigned ii = 1; ii < cnt_loci; ii++)
	{
		for (unsigned jj = 0; jj < matrix_size; jj++)
		{
			//
			opt_curr[jj]	= -1;
			path[ii][jj]	= -1;
			for (unsigned kk = 0; kk < matrix_size; kk++)
			{
				//
				//	penalise if state change (jj != k)
				//	matching = +1
				//
				//int x = opt[ii-1][kk] + mat[ii][jj] - penalty * (jj != kk);
				//if ( x > opt[ii][jj] )
				int x = opt_prev[kk] + mat[ii][jj] - penalty_matrix[jj][kk];
				if ( x > opt_curr[jj] )
				{
					opt_curr[jj]	= x;
					path[ii][jj]	= kk;
				}
			}
		}
		opt_curr.swap(opt_prev);
	}


	// 
	// find best opt at last column to start back tracing
	// p = start of optimal path
	//
	int max_opt_end = opt_prev[0];
	int p = 0;
	for (unsigned jj = 0; jj < matrix_size; jj++)
	{
		if ( opt_prev[jj] > max_opt_end )
		{
			max_opt_end = opt_prev[jj];
			p = jj;
		}
	}

	//
	//	walk backwards along optimum path
	//
	vector<unsigned>			step(cnt_loci + 1, 0);
	step[cnt_loci] = p;
	unsigned ii = cnt_loci;
	while (ii > 0)
	{
		p = path[--ii][p];
		step[ii] = p;
	} 


	//	
	//	Save bounds of intervals from a single founder (i.e. count states + 1)
	// 
	boundaries.clear();
	founder_ids.clear();
	boundaries.push_back(0);
	founder_ids.push_back(step[0]);
	for (unsigned ii=1; ii <= cnt_loci; ii++)
	{
		assert(step[ii] < matrix_size);
		if ( step[ii] != step[ii-1] )	
		{
			boundaries.push_back(ii);
			founder_ids.push_back(step[ii]);
		}
	}
	boundaries.push_back(cnt_loci);
}



//!_______________________________________________________________________________________
//!     
//!		reconstruct_mosaic_haploid
//! 
//! 		Use dynamic programming to determine which founders contributed to which
//!     		intervals of the genome
//!  		Assumes homozygosity
//!     
//!		\author lg (4/11/2011)
//!     
//!_______________________________________________________________________________________
void reconstruct_mosaic_haploid(prog_options&					opt, 
								const vector<string>& 			strain_names, 
								const string&					contig_name, 
								const t_seq_calls_per_contig& 	seq_calls, 
								const t_founder_contig&			founder_contig,
								vector<t_founder_interval>& 	founder_intervals,
								unsigned						penalty,
								unsigned						homozygous_match,
								unsigned						heterozygous_match
								)

{
	opt.verbose(2) << "Reconstruct haploid mosaic\n";
	unsigned cnt_founders 		= strain_names.size();
	unsigned matrix_size 		= cnt_founders;
	unsigned cnt_loci   		= seq_calls.pos.size();

	// 
	//  initialise matrix
	// 
	multi_array<char,2> mat	(extents[cnt_loci][matrix_size]);
	const_multi_array_ref<char, 2>	founder_alleles = founder_contig.get_diallellic_snps();
	assert(founder_alleles[0].size() == matrix_size);
	assert(seq_calls.alleles.size() == cnt_loci);
	for (unsigned ii = 0; ii < cnt_loci; ii++)
	{
		unsigned index_into_founder_contig = seq_calls.indices_into_founder_contig[ii];
		assert(index_into_founder_contig < founder_alleles.size());
		for (unsigned jj = 0; jj < matrix_size; jj++)
			mat[ii][jj] = (seq_calls.alleles[ii] == founder_alleles[index_into_founder_contig][jj]) * homozygous_match;
	}


	//
	//	Initialise penalty matrix
	//
	multi_array<unsigned,2> 	penalty_matrix(extents[matrix_size][matrix_size]);
	std::fill_n(penalty_matrix.data(), penalty_matrix.num_elements(), penalty);
	for (unsigned ii = 0; ii < matrix_size; ii++)
		penalty_matrix[ii][ii] = 0;



	//
	//	find best path through matrix 
	// 
	vector<unsigned>			boundaries;
	vector<unsigned>			founder_ids;
	find_optimal_path_via_dp(	mat,
								cnt_loci,
								matrix_size,
								boundaries,
								founder_ids, 
								penalty_matrix);
	

	//
	//	return results
	//

	assert(cnt_loci == seq_calls.indices_into_founder_contig.size());
	assert(cnt_loci == mat.size());
	assert(cnt_founders = mat[0].size());
	assert(boundaries.size() >= founder_ids.size());

	for (unsigned ii=0; ii < founder_ids.size(); ii++)
	{
		vector<unsigned>		mismatch_loci;
		unsigned founder_id	= founder_ids[ii];
		unsigned beg		= boundaries[ii];
		unsigned back 		= boundaries[ii + 1] - 1;
		assert(beg < back);
		assert(back < cnt_loci);
		assert(founder_id < cnt_founders);

		// save mismatches
		for (unsigned jj = beg; jj <= back; jj++)
			if (!mat[jj][founder_id])
				mismatch_loci.push_back(seq_calls.pos[jj]);

		auto founder_contig_beg_index  = seq_calls.indices_into_founder_contig[beg];
		auto founder_contig_back_index = seq_calls.indices_into_founder_contig[back];
		assert(founder_contig_beg_index  < founder_contig.loci.size()) ;
		assert(founder_contig_back_index < founder_contig.loci.size());
		founder_intervals.emplace_back( founder_contig_beg_index ,
										founder_contig_back_index,
										founder_contig.loci[founder_contig_beg_index ],
										founder_contig.loci[founder_contig_back_index],
										founder_id,
										mismatch_loci.size(),
										back - beg + 1);
		founder_intervals.back().mismatch_loci.swap(mismatch_loci);
	}
	// 
	// first and last intervals begin and end at first and last founder SNP position and not
	// at the first and last positions of the sparse genotypes in the strains
	// 
	// ??????????
	// 
	// 	Why does richard's code start from 1?
	// 
	// ??????????

	auto& first = founder_intervals[0];
		first.founder_contig_beg_index	= 0;
		first.beg_locus 				= founder_contig.loci[0];
	auto& last = founder_intervals.back();
		last.founder_contig_back_index	= founder_contig.cnt_sites - 1;
		last.back_locus					= founder_contig.loci[founder_contig.loci.size() - 1];
}


//!_______________________________________________________________________________________
//!     
//!		reconstruct_mosaic_diploid
//! 
//! 		Use dynamic programming to determine which founders contributed to which
//!     		intervals of the genome
//!  		allows heterozygosity
//!     
//!		\author lg (4/11/2011)
//!     
//!_______________________________________________________________________________________
void reconstruct_mosaic_diploid(prog_options&					opt, 
								const vector<string>&			strain_names, 
								const string&					contig_name, 
								const t_seq_calls_per_contig& 	seq_calls, 
								const t_founder_contig&			founder_contig,
								vector<t_founder_interval>& 	founder_intervals,
								unsigned						penalty,
								unsigned						homozygous_match,
								unsigned						heterozygous_match
								)
{
	opt.verbose(2) << "Reconstruct diploid mosaic\n";
	unsigned cnt_founders 		= strain_names.size();
	unsigned matrix_size 		= cnt_founders * (cnt_founders + 1) / 2;
	unsigned cnt_loci   		= seq_calls.pos.size();

	// 
	//  initialise matrix
	// 
	multi_array<char,2> 			mat			(extents[cnt_loci][matrix_size]);
	const_multi_array_ref<char, 2>	founder_alleles = founder_contig.get_diallellic_snps();
	assert(founder_alleles[0].size() == cnt_founders);
	assert(seq_calls.alleles.size() == cnt_loci);

	//  Heterozygous state
	enum {eEITHER, e1st, e2nd, eHET_STATES};
	// 	MATCH / MISMATCH
	multi_array<char,2> matches	(extents[cnt_loci][matrix_size]);
	enum {eUNINFORMATIVE, eMATCH_OTHER_FOUNDER, eMATCH, eMISMATCH, eMATCH_CATEGORIES};


	// 
	// 	Initialise heterozygous state transitions (i.e. match which strand)
	// 
	multi_array<char,2> het_transitions(extents[eHET_STATES][eHET_STATES]);
	for (unsigned ii = 0; ii < eHET_STATES; ++ii)
		for (unsigned jj = 0; jj < eHET_STATES; ++jj)
			if (ii == eEITHER || jj == eEITHER)
				het_transitions[ii][jj] = eUNINFORMATIVE;
			else if (ii ==jj)
				het_transitions[ii][jj] = eMATCH;
			else 
				het_transitions[ii][jj] = eMATCH_OTHER_FOUNDER;


	vector<char> curr_het_state(matrix_size, eEITHER);


	for (unsigned ii = 0; ii < cnt_loci; ii++)
	{
		unsigned index_into_founder_contig = seq_calls.indices_into_founder_contig[ii];
		assert(index_into_founder_contig < founder_alleles.size());

		unsigned hh = 0;

		//	Homozygous first
		for (unsigned jj = 0; jj < cnt_founders; jj++, hh++)
		{
			if (seq_calls.alleles[ii] == founder_alleles[index_into_founder_contig][jj])
			{
				matches[ii][hh]	= eMATCH;
				mat[ii][hh]		= homozygous_match;
			}
			else
			{
				matches[ii][hh]	= eMISMATCH;
				mat[ii][hh] = 0;
			}
		}

		// 	Heterozygous second
		//  Weights if either is same as founder, otherwise 0
		for (unsigned jj = 0; jj < cnt_founders; jj++)
			for (unsigned kk = 0; kk < jj; kk++, hh++)
			{
				// default match one or another
				mat[ii][hh]				 = heterozygous_match;
				if (seq_calls.alleles[ii] == founder_alleles[index_into_founder_contig][jj])
				{
					// 	
					// 	Genotype only matches 1st founder
					// 
					if(seq_calls.alleles[ii] != founder_alleles[index_into_founder_contig][kk])
					{
						matches[ii][hh]		= het_transitions[curr_het_state[hh]][e1st];
						curr_het_state[hh]	= e1st;
					}
					// 	
					//  Genotype matches both founders
					//  	"uninformative" genotype, don't change current het state
					// 
					else
						matches[ii][hh] 	= eUNINFORMATIVE;
				}
				// 	
				// 	Genotype only matches 2nd founder
				// 
				else if (seq_calls.alleles[ii] == founder_alleles[index_into_founder_contig][kk])
				{
						matches[ii][hh] 	= het_transitions[curr_het_state[hh]][e2nd];
						curr_het_state[hh]	= e2nd;
				}
				// 	
				// 	Genotype only matches neither founder: not heterozygous match
				// 
				else
				{
						mat[ii][hh]		= 0;
						matches[ii][hh] = eMISMATCH;
				}
			}
		assert(hh == matrix_size);
	}

	//
	//  encode pair of founders into a single id corresponding to
	//  Homozygous founders +
	//  Heterozygous in a diagonal matrix
	vector<tuple<unsigned, unsigned>> map_founder_pair_ids;
	{
		unsigned hh = 0;
		//	Homozygous first
		for (unsigned jj = 0; jj < cnt_founders; jj++)
			map_founder_pair_ids.push_back(std::make_tuple(jj, jj));
		// 	Heterozygous second
		for (unsigned jj = 0; jj < cnt_founders; jj++)
			for (unsigned kk = 0; kk < jj; kk++, hh++)
					map_founder_pair_ids.push_back(std::make_tuple(jj, kk));
	}


	//
	//	Initialise penalty matrix
	//
	multi_array<unsigned,2> 	penalty_matrix(extents[matrix_size][matrix_size]);
	for (unsigned ii=0; ii < map_founder_pair_ids.size(); ii++)
	{
		unsigned founder_id_a1, founder_id_a2;
		tie(founder_id_a1, founder_id_a2)	= map_founder_pair_ids[ii];
		for (unsigned jj=0; jj < map_founder_pair_ids.size(); jj++)
		{
			if (ii == jj)
			{
				penalty_matrix[ii][jj] = 0;
				continue;
			}
			unsigned founder_id_b1, founder_id_b2;
			tie(founder_id_b1, founder_id_b2)	= map_founder_pair_ids[jj];
			unsigned p = 0;
			if (founder_id_b1 != founder_id_a1)
				p += penalty;
			if (founder_id_b2 != founder_id_a2)
				p += penalty;
			penalty_matrix[ii][jj] = p;
		}
	}


	//
	//	find best path through matrix 
	// 
	vector<unsigned>				boundaries;
	vector<unsigned>				founder_pair_ids;
	find_optimal_path_via_dp(	mat,
								cnt_loci,
								matrix_size,
								boundaries,
								founder_pair_ids,
								penalty_matrix);



	//
	//	return results
	//
	assert(boundaries.size() == founder_pair_ids.size() + 1);
	for (unsigned ii=0; ii < founder_pair_ids.size(); ii++)
	{
		unsigned founder_pair_id		= founder_pair_ids[ii];
		unsigned founder_id1;
		unsigned founder_id2;
		assert(founder_pair_id < map_founder_pair_ids.size());
		tie(founder_id1, founder_id2)	= map_founder_pair_ids[founder_pair_id];
		unsigned beg					= boundaries[ii];
		unsigned end  					= boundaries[ii + 1];
		unsigned back 					= end - 1;

		assert(back < seq_calls.indices_into_founder_contig.size());
		assert(beg < seq_calls.indices_into_founder_contig.size());
		assert(beg < mat.size());
		assert(end <= mat.size());
		assert(founder_pair_id < mat[0].size());


		unsigned cnt_snps	= end - beg;

		// {eUNINFORMATIVE, eMATCH_OTHER_FOUNDER, eMATCH, eMISMATCH, eMATCH_CATEGORIES}
		vector<unsigned> counts_per_match_category(eMATCH_CATEGORIES, 0);
		for (unsigned jj = beg; jj < end; jj++)
			counts_per_match_category[matches[jj][founder_pair_id]] ++;
		// save errors
		vector<unsigned>		mismatch_loci;
		for (unsigned jj = beg; jj < end; jj++)
			if (matches[jj][founder_pair_id] == eMISMATCH)
				mismatch_loci.push_back(seq_calls.pos[jj]);

		unsigned cnt_mismatches 	= counts_per_match_category[eMISMATCH];
		unsigned cnt_transitions	= counts_per_match_category[eMATCH_OTHER_FOUNDER];
		unsigned cnt_informative	= counts_per_match_category[eMATCH_OTHER_FOUNDER] + counts_per_match_category[eMATCH];
		unsigned cnt_uninformative	= counts_per_match_category[eUNINFORMATIVE];

		founder_intervals.emplace_back(	seq_calls.indices_into_founder_contig[beg], 
										seq_calls.indices_into_founder_contig[back],
										seq_calls.pos[beg],
										seq_calls.pos[back],
										founder_id1,
										founder_id2,
										cnt_mismatches,
										cnt_transitions,
										cnt_informative,
										cnt_uninformative,
										cnt_snps);
		founder_intervals.back().mismatch_loci.swap(mismatch_loci);
	}
	// 
	// first and last intervals begin and end at first and last founder SNP position and not
	// at the first and last positions of the sparse genotypes in the strains
	// 

	auto& first = founder_intervals[0];
		first.founder_contig_beg_index	= 0;
		first.beg_locus 				= founder_contig.loci[0];
	auto& last = founder_intervals.back();
		last.founder_contig_back_index	= founder_contig.cnt_sites - 1;
		last.back_locus					= founder_contig.loci[founder_contig.loci.size() - 1];
}

//!_______________________________________________________________________________________
//!     
//!		open_output_path
//!     
//!		\author lg (4/27/2011)
//!     
//!     \param output_file_name 
//!     \param output_stream 
//!_______________________________________________________________________________________
void open_output_stream(const string& output_file_name, ofstream& output_stream)
{
	path output_path(output_file_name);
	auto parent_path = output_path.parent_path();
	if (!exists(parent_path))
		create_directories(parent_path);
	else if (!is_directory(parent_path))
		throw_logic_error(format("Cannot create %s: %s is not a directory.\n")
								% output_file_name % parent_path.string());
	output_stream.open(output_file_name.c_str());
	if (!output_stream)
		throw std::logic_error("Can't open " + output_file_name);
}

//!_______________________________________________________________________________________
//!     
//!		write_errors_per_contig
//!     
//!		\author lg (4/27/2011)
//!     
//!     \param output_file_name 
//!     \param mismatch_loci 
//!_______________________________________________________________________________________
void write_errors_per_contig(const string& output_file_name, vector<unsigned>& mismatch_loci)
{
	ofstream output_stream;
	open_output_stream(output_file_name, output_stream);
	for_each(mismatch_loci, [&](unsigned ii) {output_stream << ii << "\n";});
}

//!_______________________________________________________________________________________
//!     
//!		write_founders_per_contig
//!			Write the founders for each segment of each contig of each seq call file
//!     
//!		\author lg (4/21/2011)
//!     
//!     \param output_file_name 
//!     \param contig_name 
//!     \param strain_names 
//!     \param founder_intervals 
//!_______________________________________________________________________________________
void write_founders_per_contig(	ostream& 							founders_output_stream,
								const string& 						call_name,
								const string& 						contig_name,
								const string& 						haploid_or_diploid,
								const vector<string>&				strain_names,
								const vector<t_founder_interval>& 	founder_intervals)
{
	unsigned	cnt_mismatches;
	unsigned	cnt_transitions;
	unsigned	cnt_informative;
	unsigned	cnt_uninformative;
	unsigned	cnt_snps;


	BOOST_FOREACH(auto& fi, founder_intervals) 
	{
		founders_output_stream			
									<< call_name
						<< "\t"		<< contig_name 
						<< "\t"		<< haploid_or_diploid 
						<< "\t"		<< 
						((fi.founder_strain_id1 == fi.founder_strain_id2) ?
							strain_names[fi.founder_strain_id1]:
							strain_names[fi.founder_strain_id1] + "/" + strain_names[fi.founder_strain_id2])
						<< "\t" 	<< fi.beg_locus
						<< "\t" 	<< fi.back_locus
						<< "\t" 	<< fi.len()
						<< "\t" 	<< fi.cnt_snps
						<< "\t" 	<< fi.cnt_mismatches 
						<< "\t" 	<< fi.cnt_transitions
						<< "\t" 	<< fi.cnt_informative
						<< "\t" 	<< fi.cnt_uninformative
						<< "\t" 	<< fi.errors_per_snp()
						<< "\t" 	<< fi.errors_per_bp() << "\t";
		if (!fi.mismatch_loci.size())
		{
			founders_output_stream	<< "NA\n";
			continue;
		}
		founders_output_stream	<< fi.mismatch_loci[0];
		for (unsigned ss = 1, ss_end = fi.mismatch_loci.size(); ss < ss_end; ++ss)
			founders_output_stream	<< "," << fi.mismatch_loci[ss];
		founders_output_stream	<< "\n";
	}
}



//!_______________________________________________________________________________________
//! 
//!		impute_genomes
//!						
//! 		There are only a limited number of unique Strain Distribution Patterns across the genome,        
//!         At each site, we can replace the actual genotypes with                                           
//!             an index into the list of unique SDPs                                                        
//!                                                                                                          
//!         Simple dumb code:                                                                                
//!             We look at each successive site:                                                             
//!                 1) Is it a newly encountered unique SDP?                                                 
//!                     -> make a note of it                                                                 
//!                     -> link it to its own index                                                          
//!                 2) Replace the SDP with its representative index                                         
//!                                                                                                          
//!         In the original Mott code, SDP indices are in order of their string sorted value                 
//!         In the Leo code, SDP indices are in order of their first appearance in genomic order             
//!                                                                                                          
//!         We can translate from one to another for comparison...                                           
//! 
//! 
//!		\author lg (17/04/2011)
//!     
//!     \param opt 
//!     \param founder_contigs 
//!     \param founder_intervals_per_contig_per_strain 
//!     \param contig_names 
//!_______________________________________________________________________________________
//
//	hash function for matrix 
// 
template<typename CONT>  struct hash_container : public std::unary_function<CONT, std::size_t> 
{
	std::size_t operator()(CONT const& c) const
	{
		std::size_t seed = 0;
		for(auto ii = c.begin(), end = c.end(); ii != end; ++ii)
			hash_combine(seed, *ii);
		return seed;
	}
};
void impute_genomes(prog_options&										opt,
					const vector<t_founder_contig>&						founder_contigs,
					const multi_array<vector<t_founder_interval>,2>&	founder_intervals_per_contig_per_strain,
					const vector<string>&								contig_names,
					const vector<string>&								call_names)
{
	unsigned cnt_contigs		= founder_contigs.size();
	unsigned cnt_called_strains	= founder_intervals_per_contig_per_strain[0].size();
	assert(cnt_contigs == founder_intervals_per_contig_per_strain.size());
	assert(call_names.size() == cnt_called_strains);

	set<string> used_file_names;

	for (unsigned cc = 0; cc < cnt_contigs; ++cc)
	{
		//
		//	open output files 
		// 
		auto imputed_filename = replace_all_copy(opt.imputed_output_file_pattern, "[CONTIG]", contig_names[cc]);
		if (used_file_names.count(imputed_filename))
			throw_logic_error(format("File name [%s] to hold imputed genotypes has been used before.") % imputed_filename);

		ofstream imputed_stream;
		open_output_stream(imputed_filename, imputed_stream);

		auto sdp_filename = replace_all_copy(opt.sdp_output_file_pattern, "[CONTIG]", contig_names[cc]);
		if (used_file_names.count(sdp_filename))
			throw_logic_error(format("File name [%s] to hold strain distribution patterns has been used before.") % sdp_filename);
		ofstream sdp_stream;
		open_output_stream(sdp_filename, sdp_stream);

		opt.log() << format("Chr%1$s Imputed indices / SDP in [%2$s / %3$s]\n") % contig_names[cc] % imputed_filename % sdp_filename;

		//
		//	local references to the founder alleles for this contig
		// 
		const t_founder_contig& founder_contig		= founder_contigs[cc];
		const_multi_array_ref<char, 2>	sdp			= founder_contig.get_sdp();
		unsigned cnt_sites							= founder_contig.cnt_sites;



		// ????????????????????????????????????????
		// 
		// why + 1?? in Richard Mott's code
		// 
		// ????????????????????????????????????????
		multi_array<char, 2> 	mat(extents[cnt_sites][cnt_called_strains]);
		vector<unsigned>		ptr(make_counting_iterator(0U), make_counting_iterator(cnt_sites));


		//
		//  Assign SDP of founders to matrix of strains using imputed haplotypes
		// 
		opt.verbose(4) << "Assign SDP of founders to matrix of strains using imputed haplotypes\n";
		assert(cnt_called_strains == mat[0].size());
		for (unsigned ss = 0; ss < cnt_called_strains; ++ss)
		{
			const vector<t_founder_interval>& founder_intervals = founder_intervals_per_contig_per_strain[cc][ss];
			for (unsigned kk = 0, kk_end = founder_intervals.size(); kk < kk_end; ++kk)
			{
				unsigned from_idx 			= founder_intervals[kk].founder_contig_beg_index;  
				unsigned to_idx				= founder_intervals[kk].founder_contig_back_index;  
				unsigned founder_strain_id	= founder_intervals[kk].founder_strain_id1;
				assert(founder_strain_id < sdp[0].size()) ;
				if (from_idx >= sdp.size())
					std::cerr << "Oops from" << from_idx << ", " << sdp.size() << "\n";
				if (to_idx >= sdp.size())
					std::cerr << "Oops to" << to_idx << ", " << sdp.size() << "\n";
				assert(from_idx < sdp.size());
				assert(to_idx < sdp.size()) ;
				assert(from_idx < mat.size()) ;
				assert(to_idx < mat.size()) ;
				for (unsigned mm = from_idx; mm <= to_idx; ++mm)
					mat[mm][ss] = sdp[mm][founder_strain_id];
			}
		}

		opt.verbose(4) << "Replace genotypes with an index into the list of unique SDPs\n";
		//
		//  There are only a limited number of unique Strain Distribution Patterns across the genome,
		// 	At each site, we can replace the actual genotypes with 
		// 		an index into the list of unique SDPs
		// 

		// hash from SDP to its index in the list of unique SDPs
		typedef const_multi_array_ref<char, 2>::value_type sub_array_type;
		typedef unordered_map<sub_array_type, unsigned, hash_container<sub_array_type>> map_array_u;
		map_array_u map_sdp_to_index;             
		vector<unsigned> sdp_index_per_site;

		// the offset into the list of sites (in "mat") where the successive unique SDPs can be found
		vector<unsigned> unique_sdp_indices;

		sdp_index_per_site.reserve(mat.size());
		map_uu unique_sdp_index_popularity;
		for (unsigned mm = 0, end = mat.size(); mm < end; ++mm)
		{
			// 
			// only change map_sdp_to_index if it wasn't there already
			//
			auto ret = map_sdp_to_index.insert(make_pair(mat[mm], map_sdp_to_index.size()));
			// in which case save the index of this newly encountered sdp


			if (ret.second)
				unique_sdp_indices.push_back(mm);

			// replace the entire SDP with just its index
			auto iter = ret.first;
			sdp_index_per_site.push_back(iter->second);
			++unique_sdp_index_popularity[iter->second];
		}
		opt.verbose(4) << "Output imputed SDP indices\n";


		//
		//	sort SDP indices by popularity
		// 
		vector<std::pair<unsigned, unsigned>> vec_popularity_sdp_index;
		BOOST_FOREACH(auto& p, unique_sdp_index_popularity) 
			vec_popularity_sdp_index.push_back({p.second, p.first});
		sort(vec_popularity_sdp_index);
		reverse(vec_popularity_sdp_index);

		map_uu sdp_index_to_by_popularity;
		for (unsigned ii = 0, ii_end = vec_popularity_sdp_index.size(); ii < ii_end; ++ii)
			sdp_index_to_by_popularity[vec_popularity_sdp_index[ii].second] = ii;

		
		// 
		// output imputed SDP indices
		// 
		assert(sdp_index_per_site.size() == founder_contigs[cc].loci.size());
		for (unsigned ss = 0, ss_end = sdp_index_per_site.size(); ss < ss_end; ++ss)
		{
			unsigned sdp_index_by_popularity = sdp_index_to_by_popularity[sdp_index_per_site[ss]];
			imputed_stream	<< contig_names[cc] << "\t" 
							<< founder_contigs[cc].loci[ss] << "\t"
							<< sdp_index_by_popularity << "\n";
		}

		// 
		// 	Output actual SDPs
		// 
		// header
		opt.verbose(4) << "Output actual SDPs\n";
		sdp_stream << "SDP";
		for (unsigned ss = 0, ss_end = vec_popularity_sdp_index.size(); ss < ss_end; ++ss)
			sdp_stream << "\t" << ss;
		sdp_stream << "\n";

		// 
		// 	Each SDP in a column with each row == the genotype at a called strain
		// 
		for (unsigned cs = 0; cs < cnt_called_strains; ++cs)
		{
			sdp_stream << call_names[cs];
			BOOST_FOREACH(auto& ss, vec_popularity_sdp_index) 
				sdp_stream << "\t" << static_cast<unsigned>(mat[unique_sdp_indices[ss.second]][cs]);
			sdp_stream << "\n";
		}



		// 
		// popularity of each sdp, and its genotype pattern
		//
		opt.summary() << "SDP index" << "\t" << "#Loci with SDP" << "\t" << "Genotype Pattern\n";
		for (unsigned ss = 0, ss_end = vec_popularity_sdp_index.size(); ss < ss_end; ++ss)
		{
			unsigned	sdp_index = unique_sdp_indices[vec_popularity_sdp_index[ss].second];
			map_uu		cnt_genotypes;
			for (unsigned cs = 0; cs < cnt_called_strains; ++cs)
				++cnt_genotypes[mat[sdp_index][cs]];
			opt.summary() << ss << "\t" << vec_popularity_sdp_index[ss].first << "\t";
			BOOST_FOREACH(auto& gg, cnt_genotypes) 
				opt.summary() << gg.first << "=" << gg.second << "; ";
			opt.summary() << "\n";
		}
	}
}


//!_______________________________________________________________________________________
//!     
//!		get_call_name
//! 		Construct name to tag a set of calls from the call file name using
//! 		a regular expression
//!     
//!		\author lg (4/21/2011)
//!     
//!     \param opt 
//!     \param call_file_to_name_regex 
//!     \param file_name 
//!     \param set_of_call_names 
//!     
//!     \return string 
//!_______________________________________________________________________________________
string get_call_name(	prog_options& 			opt, 
						const boost::regex&		call_file_to_name_regex, 
						const string&			file_name, 
						set<string>&			set_of_call_names)
{
	//
	//  Construct unique name to identify data from this call file
	// 
	const string& call_name =  	boost::regex_replace(file_name, call_file_to_name_regex, "$1", boost::match_default | boost::format_all);
	if (!call_name.length())
		throw_logic_error(format("Unable to construct a name for this sequence call file [%s] "
								 "using the regular expression --call_file_to_name_regex '%s'. "
								 "Are you using a valid regular expression? Have you forgotten to "
								 "add brackets around the name?")
									% file_name % opt.call_file_to_name_regex);
	if (set_of_call_names.count(call_name))
		throw_logic_error(format("The name [%s] does not uniquely identify this sequence call file [%s] "
								 "using the regular expression --call_file_to_name_regex '%s'. "
								 "Are you using a valid regular expression?")
									% call_name % file_name % opt.call_file_to_name_regex);
	set_of_call_names.insert(call_name);
	return call_name;
}


//!_______________________________________________________________________________________
//!     
//!		get_output_file_name_per_contig
//!			
//!     
//!		\author lg (4/21/2011)
//!     
//!     \param opt 
//!     \param file_name 
//!     \param contig_name 
//!     \param set_of_call_names 
//!     
//!     \return string 
//!_______________________________________________________________________________________
string get_output_file_name_per_contig(	prog_options& 	opt, 
										const string&	output_file_pattern, 
										const string&	pattern_name, 
										const string&	file_name, 
										const string&	call_name, 
										const string&	contig_name, 
										set<string>&	set_of_output_names)
{
	//
	//  Construct unique name to identify data from this call file
	// 

	std::string output_file_name = get_file_name_substitutions(output_file_pattern,
																		file_name,
																		{{"[CALLNAME]",	call_name},
																		{"[CONTIG]", 	contig_name}});
	if (!output_file_name.length())
		throw_logic_error(format("Unable to construct an output file name for "
								 "this sequence call file [%s] using "
								 "%s '%s'. Are you using a valid pattern?")
									% file_name % pattern_name % output_file_pattern);
	if (set_of_output_names.count(output_file_name))
		throw_logic_error(format("The output file name [%s] is not unique and has been used already for data from [%s] "
								 "specified by '%s '%s'. \n"
								 "Are you using a valid pattern containing '[FILENAME]' or '[CALLNAME]'?")
									% output_file_name % file_name % pattern_name % output_file_pattern);

	set_of_output_names.insert(output_file_name);
	return output_file_name;
}


//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
// 
//      do_main
// 
//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888 
void do_main(prog_options& opt)
{

	// 
	//  Read regions to exclude
	// 
	t_intervals_per_contig excluded_intervals_by_contig;
	BOOST_FOREACH(const string& file_name, opt.exclude_regions_gff_filenames) 
	{
		t_intervals_per_contig curr_excluded_intervals_by_contig;
		read_intervals(	opt, file_name, curr_excluded_intervals_by_contig, 
								false, //zero_based_coordinates, 
								intervals_per_contig::eGFF_FORMAT);
		excluded_intervals_by_contig.combine_with(curr_excluded_intervals_by_contig);
	}
	BOOST_FOREACH(const string& file_name, opt.exclude_regions_filenames) 
	{
		t_intervals_per_contig curr_excluded_intervals_by_contig;
		read_intervals(	opt, file_name, curr_excluded_intervals_by_contig, 
								false, //zero_based_coordinates, 
								intervals_per_contig::eRAW_FORMAT);
		excluded_intervals_by_contig.combine_with(curr_excluded_intervals_by_contig);
	}

	//
	//  Read founder alleles
	// 
	vector<t_founder_contig>	founder_contigs;
	vector<string>				contig_names;
	vector<string>				strain_names;
	BOOST_FOREACH(const string& file_name, opt.allele_filenames) 
		read_allele_data(opt, file_name, founder_contigs, contig_names, strain_names);



	//
	//	Map genomic positions to alleles
	// 
	vector<map_uu>	map_pos_to_allele_index_per_contig(contig_names.size(), map_uu());



	//
	//	map positions to founder alleles index, excluding specified regions
	// 		if no specified regions, will just map (i.e. so quite fast)
	// 	Only founder alleles which are diallelic SNPs (ok == true)
	//
	for (unsigned cc = 0; cc < contig_names.size(); ++cc)
	{
		
		// 
		// 	No matching contig among exclusion intervals.
		//  Create map directly
		//
		ci_map_str_intervals ff = excluded_intervals_by_contig.data.find(contig_names[cc]);
		map_uu&	contig_map = map_pos_to_allele_index_per_contig[cc];
		t_founder_contig& contig_alleles = founder_contigs[cc];
		const string& contig_name = contig_names[cc];

		if (ff == excluded_intervals_by_contig.data.end())
		{
			opt.verbose(2) << contig_name << " has no excluded intervals.\n";
			for (unsigned jj = 0; jj < contig_alleles.loci.size(); ++jj)
			{
				if (!contig_alleles.ok[jj])
					continue;
				contig_map[contig_alleles.loci[jj]] = jj;
			}
		}

		// 
		// 	Check each locus is diallelic SNP and not in excluded regions before mapping
		//
		else
		{
			t_interval_search_cache cache_for_search;
			const t_intervals& exclusion = ff->second;
			for (unsigned jj = 0; jj < contig_alleles.loci.size(); ++jj)
			{
				if (!contig_alleles.ok[jj])
					continue;
				if (!exclusion.contains(contig_alleles.loci[jj], cache_for_search))
					contig_map[contig_alleles.loci[jj]] = jj;
			}
		}
		opt.both(3) << format("    Chromosome %s: %12s mapped alleles falling outside excluded intervals.\n") 
						% contig_name 
						% contig_map.size();
	}


	// 
	//  Read sequence calls
	// 
	opt.both() << format("%d files with sequence calls.\n") % opt.call_filenames.size();
	vector<t_seq_calls_per_contig>		seq_calls_per_contig;

	//
	//	returns founder intervals per contig per seq call
	// 


	boost::regex call_file_to_name_regex(opt.call_file_to_name_regex);

	multi_array<vector<t_founder_interval>,2> 	haploid_founder_intervals_per_contig_per_strain (extents[contig_names.size()][opt.call_filenames.size()]);
	multi_array<vector<t_founder_interval>,2> 	diploid_founder_intervals_per_contig_per_strain (extents[contig_names.size()][opt.call_filenames.size()]);
	set<string> set_of_call_names;
	set<string> set_of_output_names;
	vector<string> call_names;
	for (unsigned ss = 0; ss < opt.call_filenames.size(); ++ss)
	{
		const string& call_file_name =  opt.call_filenames[ss];

		//
		//  Construct unique name to identify data from this call file
		// 
		string call_name =  get_call_name(	opt, 
											call_file_to_name_regex, 
											call_file_name, 
											set_of_call_names);
		call_names.push_back(call_name);

		

		t_portable_timer	timer;
		read_seq_calls(	opt, 
						call_file_name, 
						founder_contigs.size(),
						map_pos_to_allele_index_per_contig,
						seq_calls_per_contig,
						contig_names);

		assert(seq_calls_per_contig.size() == contig_names.size());
		assert(founder_contigs.size() == contig_names.size());

		ofstream founders_output_stream;
		open_output_stream(opt.founders_output_file, founders_output_stream);
		// 
		// write header
		//
		founders_output_stream		<< "call_name"
							"\t" 	"chr" 
							"\t" 	"haploid_diploid" 
							"\t" 	"founders"
							"\t" 	"from.bp"
							"\t" 	"to.bp"
							"\t" 	"len"
							"\t" 	"snps"
							"\t" 	"mismatches"
							"\t" 	"transitions"
							"\t" 	"informative"
							"\t" 	"uninformative"
							"\t" 	"mismatches_per_SNP"
							"\t" 	"mismatches_per_bp"
							"\t" 	"mismatch_loci"
							"\n";

		// 
		//  Reconstruct mosaic
		// 
		for (unsigned cc = 0; cc < contig_names.size(); ++cc)
		{
			if (opt.haploid)
			{
				vector<unsigned> mismatch_loci;
				reconstruct_mosaic_haploid(	opt, strain_names, contig_names[cc], 
											seq_calls_per_contig[cc], founder_contigs[cc], 
											haploid_founder_intervals_per_contig_per_strain[cc][ss], 
											opt.penalty, opt.homozygous_match, opt.heterozygous_match);
				write_founders_per_contig(	founders_output_stream,
											call_name,
											contig_names[cc],
											"haploid",
											strain_names,
											haploid_founder_intervals_per_contig_per_strain[cc][ss]);
			}

			if (opt.diploid)
			{
				vector<unsigned> mismatch_loci;
				reconstruct_mosaic_diploid(	opt, strain_names, contig_names[cc], 
											seq_calls_per_contig[cc], founder_contigs[cc], 
											diploid_founder_intervals_per_contig_per_strain[cc][ss], 
											opt.penalty, opt.homozygous_match, opt.heterozygous_match);
				write_founders_per_contig(	founders_output_stream,
											call_name,
											contig_names[cc],
											"diploid",
											strain_names,
											diploid_founder_intervals_per_contig_per_strain[cc][ss]);
			}


		}
		string time_taken = timer.pretty_time(timer.reset());
		opt.verbose() << format("    Parsing sparse SNP calls from %s in %s\n") % path(call_file_name).filename().native() % time_taken;
		opt.log() << format("Sparse SNP calls from %s in %s\n") % call_file_name % time_taken;
	}

	impute_genomes(	opt,
					founder_contigs,
					haploid_founder_intervals_per_contig_per_strain,
					contig_names,
					call_names);

	opt.verbose() << "Done\n";
}



